function [Z,W,eta,invlambda,zcidx,newK,etime,fval] = mcmc(Z, W, eta, invlambda, X, y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, C, ell, poisstrunc, algtype, loopi)
[Nall,K] = size(Z);
[D,tK] = size(W);
Ntrain = numel(invlambda);
if K ~= tK || K ~= numel(eta) || Nall~=size(X,2) || D~=size(X,1) || Ntrain~=numel(y)
    error('Please check data dimension.');
end

if ~any(algtype == 1:2)
    error('cfg:invalidparam', ...
        ['Unknown ''algtype'': %d\n', ...
         '    1 for ''semi-collapsed Gibbs sampling'';\n', ...
         '    2 for ''semi-collapsed EM-style''.'], algtype);
end

etime = 0;
ellp = 0.5*C*ell;
yp = 0.5*C*y;
for e = 1:1
    % sampling invlambda
    tStart = tic;
    tmu = 1./abs(ellp-yp.*(Z(1:Ntrain,:)*eta));
    indinf = isinf(tmu);
    if any(indinf(:))
        tmu(indinf) = max(tmu(~indinf)).^2; % avoid infinite mu
    end
    if algtype == 2
        invlambda = tmu;
    else
        invlambda = invnrnd(tmu, 1);
    end
    tElapsed = toc(tStart);
    fprintf('%.4f[%d]: %s (%.2fs) | invlambda\n', C, loopi, ...
        num2str([mean(invlambda(:)),std(invlambda(:),1)], ['%.2f',char(177),'%.2f']), tElapsed);
    etime = etime + tElapsed;
    
    % sampling W
    tStart = tic;
    tmp = sqrt(invsigmaxsqr);
    tmp = Z.*tmp(:,ones(K,1));
    invU = tmp'*tmp;
    invU(1:K+1:K*K) = invU(1:K+1:K*K) + invsigmawsqr;
    R = choll(invU);
    for d = 1:D
        ud = invsigmaxsqr.*X(d,:)';
        ud = R\(R'\(sum(Z.*ud(:,ones(1,K)))'));
        if algtype == 2
            W(d,:) = ud';
        else
            W(d,:) = (ud + R\randn(K, 1))';
        end
    end
    tElapsed = toc(tStart);
    fval = fobj(Z, W, eta, X, y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, C, ell);
    fprintf('%.4f[%d]: %s (%.2fs) | W\n', C, loopi, ...
        num2str([mean(W(:)),std(W(:),1),fval], ['%.2f',char(177),'%.2f',', %.4f']), tElapsed);
    etime = etime + tElapsed;
    
    % sampling Z
    tStart = tic;
    newks = zeros(1,poisstrunc+1);
    m = sum(Z);
    for n = randperm(Nall)
        m = m - Z(n,:);
        WZnt = W*Z(n,:)';
%         for k = 1:K
%             oZik = Z(n,k);
%             p1bp0 = m(k)/(Nall-m(k));
%             p1bp0 = p1bp0 * exp(-invsigmaxsqr(n)*(W(:,k)'*(WZnt-X(:,n)+(0.5-oZik)*W(:,k))));
%             if n <= Ntrain
%                 Z(n,k) = 0;
%                 tmp = ellp-yp(n)*Z(n,:)*eta+1./invlambda(n);
%                 p1bp0 = p1bp0 * exp(-0.5*((tmp-yp(n)*eta(k))^2-tmp^2)*invlambda(n));
%             end
%             rval = rand;
%             if algtype == 2
%                 Z(n,k) = p1bp0 > 1;
%             else
%                 Z(n,k) = rval*(1+p1bp0) > 1;
%             end
% %             fprintf('%.3f\t%.3f\t%d\n', p1bp0, rval, Z(n,k));
%             WZnt = WZnt+(Z(n,k)-oZik)*W(:,k);
%         end
        zrow = Z(n,:);
        if n > Ntrain
            zrow_sampler(zrow, m./(Nall-m), W, -X(:,n), ones(D,1), WZnt, invsigmaxsqr(n), algtype);
        else
            zrow_sampler(zrow, m./(Nall-m), W, -X(:,n), ones(D,1), WZnt, invsigmaxsqr(n), algtype, ...
                yp(n)*eta, -ellp-1./invlambda(n), yp(n)*zrow*eta, invlambda(n));
        end
        Z(n,:) = zrow;
        m = m + Z(n,:);
        
        % semi-collapsed sampling
        [detsA, suminvA] = scgibbshelper(invsigmawsqr, invsigmaxsqr(n), poisstrunc);
        loglhdk = 0.5*D*((1:poisstrunc)*log(invsigmawsqr) - log(detsA)) + (0.5*invsigmaxsqr(n)^2*sum((X(:,n)-WZnt).^2)*ones(1,poisstrunc)).*suminvA;
        if n <= Ntrain
            [detsB, suminvB] = scgibbshelper(invsigmaetasqr, 0.25*C^2*invlambda(n), poisstrunc);
            tc = yp(n)*(1+(ellp-yp(n)*(Z(n,:)*eta))*invlambda(n));
            loglhdk = loglhdk + 0.5*((1:poisstrunc)*log(invsigmaetasqr) - log(detsB) + (tc^2*ones(1,poisstrunc)).*suminvB);
        end
        logpk = [0, cumsum(log(alphav/N)-log(1:poisstrunc))] + [0, loglhdk];
        if algtype == 2
            [~, newk] = max(logpk);
            newk = newk - 1;
        else % algtype == 1
            newk = logmnrnd(logpk) - 1;
        end
        if newk > 0
            Z(n,K+1:K+newk) = 1;
            m(K+1:K+newk) = 1;
            tmu = invsigmaxsqr(n)*(X(:,n)-WZnt)*(suminvA(:,newk)./newk)';
            tmu = tmu(:, ones(1,newk));
            tSigma = scgibbshelperinv(detsA, suminvA, newk);
            if algtype == 2
                W(:,K+1:K+newk) = tmu;
            else % algtype == 1
                W(:,K+1:K+newk) = mvnrnd(tmu, tSigma);
            end
            if n <= Ntrain
                tmu = tc.*suminvB(:,newk)./newk;
                tmu = tmu(:, ones(1,newk));
                tSigma = scgibbshelperinv(detsB, suminvB, newk);
                if algtype == 2
                    eta(K+1:K+newk,:) = tmu;
                else
                    eta(K+1:K+newk,:) = mvnrnd(tmu, tSigma);
                end
            else
                if algtype == 2
                    eta(K+1:K+newk,:) = 0;
                else
                    eta(K+1:K+newk,:) = randn(newk, 1)./sqrt(invsigmaetasqr);
                end
            end
            K = K + newk;
        end
        newks(newk+1) = newks(newk+1) + 1;
    end
    newK = sum(newks(2:end).*(1:poisstrunc));
    zcidx = sum(Z(:,1:K-newK))==0;
    Z(:,zcidx) = [];
    W(:,zcidx) = [];
    eta(zcidx,:) = [];
    K = size(Z,2);
    tElapsed = toc(tStart);
    fval = fobj(Z, W, eta, X, y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, C, ell);
    nfeapu = sum(Z,2);
    maxnewkid = find(newks, 1, 'last');
    fprintf('%.4f[%d]: %d(%d), %s, [%s][%s], %.4f (%.2fs) | Z\n', C, loopi, ...
        K, maxnewkid-1, num2str([mean(nfeapu),std(nfeapu(:),1)], ['%.2f',char(177),'%.2f']), ...
        num2str(newks(2:maxnewkid)./(sum(newks)-newks(1)), '%.2f '), ...
        num2str(newks(1:maxnewkid)./(sum(newks)-cumsum([0,newks(1:maxnewkid-1)])), '%.2f '), ...
        fval, tElapsed);
    etime = etime + tElapsed;
    
    % sampling eta
    tStart = tic;
    tmp = 0.5*C*sqrt(invlambda);
    tmp = Z(1:Ntrain,:).*tmp(:,ones(K,1));
    invB = tmp'*tmp;
    invB(1:K+1:K*K) = invB(1:K+1:K*K) + invsigmaetasqr;
    R = choll(invB);
    b = yp.*(1+ellp.*invlambda);
    b = R\(R'\(sum(Z(1:Ntrain,:).*b(:,ones(1,K)))'));
    if algtype == 2
        eta = b;
    else
        eta = b + R\randn(K, 1);
    end
    tElapsed = toc(tStart);
    fval = fobj(Z, W, eta, X, y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, C, ell);
    fprintf('%.4f[%d]: %s (%.2fs) | eta\n', C, loopi, ...
        num2str([mean(eta(:)),std(eta(:),1),fval], ['%.2f',char(177),'%.2f',', %.4f']), tElapsed);
    etime = etime + tElapsed;
end

end